C  MODPATH release: Version 4.00 (V4, Release 1, 2-2000)
C    Changes to work with MODFLOW-2000 -- budget routines were changed and
C    moved to file BUDGETRD.FOR
C
C MODPATH Version 3.00 (V3, Release 2, 5-99)
C Changes:
C   Bug Fixes: 8-97
C     1. Subroutines RDBUDG was modified to correctly read the MODFLOW
C        head output file when head is saved as a single 2D cross section
C        using the MODFLOW-96 option, XSECTION.
C
C   Bug Fixes: 5-25-99
C     1. Fixed error in subroutine GENPAR to correctly dimension array YMAX
C
C Previous release: MODPATH Version 3.00 (V3, Release 1, 9-94)
C***** SUBROUTINES *****
C     GETLAY
C     OPNFIL
C     CHOP
C     UPCASE
C     UREADI
C     ZERO
C     BALNCE
C     CELDAT
C     FLIPQ
C     GETPS
C     GETREQ
C     GETSTP
C     STPSIZ
C     GENPAR
C     GETABT
C     ND2IJK
C     RELEAS
C     ZCALC
C***********************
 
C***** SUBROUTINE *****
      SUBROUTINE GETLAY(BUFF,KSTP,KPER,PERTIM,TOTIM,TXT,NC,NR,K,NCOL,
     1            NROW,IUHED,IO,IEND,IOP)
      DIMENSION BUFF(NCOL,NROW)
      CHARACTER*(*) TXT
      CHARACTER*132 LINE
      CHARACTER*20 FMT
      CHARACTER*16 STRING,TEXT
 
      TEXT= TXT
      CALL UPCASE(TEXT)
      CALL CHOP(TEXT,LTEXT)
      IF(LTEXT.EQ.0) LTEXT=1
 
C... IUHED < 0 MEANS READ AN UNFORMATTED HEAD FILE
      IF(IUHED.LT.0) THEN
      IU= -IUHED
C      READ(IU,END=10,ERR=30) KSTP,KPER,PERTIM,TOTIM,STRING,NC,NR,K
      READ(IU,ERR=10) KSTP,KPER,PERTIM,TOTIM,STRING,NC,NR,K
      N= -1
      IF(TEXT.NE.' ') N= INDEX(STRING,TEXT(1:LTEXT))
      IF(N.EQ.0) GO TO 20
      READ(IU) ((BUFF(J,I),J=1,NCOL),I=1,NROW)
 
C... OTHERWISE READ A FORMATTED FILE
      ELSE
      IU=IUHED
 
C... IF THIS IS A FORMATTED HEAD FILE FROM MODFLOW, READ THE HEADER
      IF(IOP.EQ.0) THEN
      READ(IU,'(A)',END=10,ERR=40) LINE
      ICOL=1
      CALL URWORD(LINE,ICOL,IWFRST,IWLAST,2,KSTP,RDUMMY,IO,IU)
      CALL URWORD(LINE,ICOL,IWFRST,IWLAST,2,KPER,RDUMMY,IO,IU)
      CALL URWORD(LINE,ICOL,IWFRST,IWLAST,3,IDUMMY,PERTIM,IO,IU)
      CALL URWORD(LINE,ICOL,IWFRST,IWLAST,3,IDUMMY,TOTIM,IO,IU)
      CALL URWORD(LINE,ICOL,IWFRST,IWLAST,1,IDUMMY,RDUMMY,IO,IU)
        IF(TEXT.NE.' ') THEN
          N= INDEX(LINE(IWFRST:IWLAST),TEXT(1:LTEXT))
          IF(N.EQ.0) GO TO 20
        END IF
      CALL URWORD(LINE,ICOL,IWFRST,IWLAST,2,NC,RDUMMY,IO,IU)
      CALL URWORD(LINE,ICOL,IWFRST,IWLAST,2,NR,RDUMMY,IO,IU)
      CALL URWORD(LINE,ICOL,IWFRST,IWLAST,2,K,RDUMMY,IO,IU)
      CALL URWORD(LINE,ICOL,IWFRST,IWLAST,1,IDUMMY,RDUMMY,IO,IU)
      FMT= LINE(IWFRST:IWLAST)
      IF(FMT.EQ.' ') FMT='(FREE)'
      IF(NC.NE.NCOL .OR. NR.NE.NROW) THEN
        WRITE(IO,*) 'INCONSISTENT DIMENSIONS IN HEAD FILE. STOP.'
		call stopfile  ! emrl
        STOP
      END IF
 
C... OR, IF THIS IS JUST A TEXT FILE THAT CONTAINS AN ARRAY, FIGURE OUT
C    WHETHER IT SHOULD BE READ FREE OR WITH FORMAT. THEN SKIP DOWN TO THE
C    LINE WHERE THE ARRAY STARTS.
      ELSE
        REWIND(IU)
        IF(TEXT.EQ.' '.OR.TEXT.EQ.'(*)'.OR.TEXT.EQ.'(FREE)') THEN
          FMT= '(FREE)'
        ELSE
          FMT=TEXT
        END IF
        IF(IOP.GT.1) THEN
          DO 1 N=1,IOP-1
          READ(IU,*,END=50)
1         CONTINUE
        END IF
      END IF
 
C... NOW READ THE ARRAY FROM THE TEXT FILE
      DO 5 I=1,NROW
        IF(FMT.EQ.'(FREE)') THEN
          READ(IU,*,END=50) (BUFF(J,I),J=1,NCOL)
        ELSE
          READ(IU,FMT,END=50) (BUFF(J,I),J=1,NCOL)
        END IF
5     CONTINUE
      END IF
c
      IEND=0
      RETURN
 
C... HANDLE ERRORS AND SPECIAL RETURN CODES
 
C...  END OF FILE REACHED
10    IEND=1
      RETURN
 
C... POSSIBLY A VALID HEADER, BUT DOES NOT HAVE "HEAD" FLAG IN TEXT STRING
20    WRITE(IO,1000) TEXT(1:LTEXT)
1000  FORMAT(1X,A,' FILE DOES NOT HAVE VALID HEADER RECORD. STOP.')
	call stopfile  ! emrl
      STOP
 
C... PROBLEM READING HEADER RECORD
30    WRITE(IO,1100) TEXT(1:LTEXT)
1100  FORMAT(' ERROR READING HEADER RECORD IN UNFORMATTED ',A,
     1' FILE. STOP.')
		call stopfile  ! emrl
      STOP
 
40    WRITE(IO,1200) TEXT(1:LTEXT)
1200  FORMAT(' ERROR READING HEADER RECORD IN FORMATTED ',A,
     1' FILE. STOP.')
		call stopfile  ! emrl
      STOP
 
50    WRITE(IO,1300) IU
1300  FORMAT(1X,'END-OF-FILE ON UNIT ',I3,'. STOP.')
	call stopfile  ! emrl
      STOP
      END
 
C***** SUBROUTINE *****
      SUBROUTINE OPNFIL (IU,FNAME,NSTAT,IOUT,IBATCH,NACT)
      CHARACTER*256 FMT,FM,ACS  ! emrl 80 to 256
      CHARACTER*(*) FNAME
	character*256 fname2,path  !emrl jig  ! emrl 80 to 256
      LOGICAL*4 EX
      INCLUDE 'openspec.inc'
	common /emrlpath/ path  ! emrl jig
C
C                     VARIABLES
C
C  IU = FORTRAN UNIT NUMBER OF FILE
C  FNAME = FILE NAME
C  NSTAT = STATUS OF FILE
C          1 => OLD FILE
C          2 => NEW FILE
C          3 => SCRATCH FILE (DELETED AUTOMATICALLY WHEN RUN ENDS)
C          4 => UNDETERMINED STATUS (MAY OR MAY NOT EXIST)
C                 IF IT DOES NOT EXIST, IT IS CREATED BY OPEN STATEMENT
C                 IF IT DOES EXIST, IT IS OPENED AS 'OLD' FILE
C
C  IOUT = UNIT NUMBER FOR OUTPUT FILE TO WRITE ERROR MESSAGE TO
C         IF NECESSARY
C
C  IBATCH = FLAG INDICATING IF THERE IS INTERACTIVE DIALOGUE AT
C           TERMINAL
C
C           0 => THERE IS INTERACTIVE DIALOGUE
C           1 => THERE IS NOT INTERACTIVE DIALOGUE (BATCH MODE)
C
C  NACT = VARIABLE DENOTING IF A FILE IS READ ONLY, WRITE ONLY, OR
C         READ AND WRITE.
C
C           1 => READ ONLY
C           2 => WRITE ONLY
C           3 => READ AND WRITE
C
C  "NACT" IS NOT USED IN THIS SUBROUTINE. ALL FILES ARE OPENED FOR
C  READING AND WRITING. OPENING "READ ONLY" AND "WRITE ONLY" FILES
C  IS MACHINE DEPENDENT. THE VARIABLE "NACT" IS PASSED TO THIS
C  ROUTINE TO MAKE IT EASY TO MODIFY THE OPEN STATEMENTS TO ALLOW
C  FOR READ AND WRITE ONLY FILES ON ANY GIVEN MACHINE. THE CALLS TO
C  THIS SUBROUTINE ARE CURRENTLY SET UP SO THAT ANY FILE THAT IS
C  WRITTEN TO IS GIVEN "NACT = 3" AND ANY FILE THAT IS READ FROM BUT
C  NOT WRITTEN TO IS GIVEN "NACT = 1". NO FILES ARE GIVEN "WRITE ONLY"
C  STATUS. "READ ONLY" STATUS IS USEFUL IF A NUMBER OF USERS WANT
C  TO SIMULTANEOUSLY SHARE INPUT FILES.
C
C  IF THE FILE DOES NOT CURRENTLY EXIST, A NEGATIVE UNIT NUMBER IS
C  USED AS A FLAG TO INDICATED THAT A NEW FILE SHOULD BE CREATED AS
C  AN UNFORMATTED (BINARY) FILE.
C
C      call setpath(path,fname)  ! emrl jig
      IBIN=0
      IF (IU.LT.0) THEN
      IU= -IU
      IBIN=1
      END IF
C
C  CHECK TO SEE IF A FILE IS OPENED TO UNIT=IOUT SO THAT ERROR
C  MESSAGES CAN BE WRITTEN.
      INQUIRE (UNIT=IOUT,OPENED=EX)
      IO=0
      IF (EX) IO=1
C
C  OPEN AN EXISTING FILE
C
      IF (NSTAT.EQ.1) THEN
10    INQUIRE (FILE=FNAME,EXIST=EX,UNFORMATTED=FM)
C  CHECK TO SEE IF FILE EXISTS
      IF (EX) GO TO 20
      IF (IBATCH.EQ.0) THEN
      WRITE (*,*) 'THE FOLLOWING FILE DOES NOT EXIST:'
      WRITE (*,'(1X,A)') FNAME
      WRITE (*,*) 'ENTER THE NAME OF AN EXISTING FILE (<CR>=QUIT):'
      READ (*,'(A)') FNAME
      IF (FNAME.EQ.' ') STOP
      GO TO 10
      ELSE
      IF (IO.EQ.1) WRITE (IOUT,*) 'FILE DOES NOT EXIST:'
      IF (IO.EQ.1) WRITE (IOUT,'(A)') FNAME
		call stopfile  ! emrl
      STOP
      END IF
20    CONTINUE
      FMT='FORMATTED'
      ACS='SEQUENTIAL'
      IF (IBIN.EQ.1.OR.FM.EQ.'YES') THEN
         FMT=FORM
         ACS=ACCESS
      END IF
      OPEN (IU,FILE=FNAME,STATUS='OLD',FORM=FMT,ACCESS=ACS,IOSTAT=IERR)
      IF (IERR.GT.0) GO TO 40
      RETURN
      END IF
C
C  OPEN A NEW FILE
C
      IF (NSTAT.EQ.2) THEN
30    INQUIRE (FILE=FNAME,EXIST=EX)
      IF (EX) THEN
      IF (IBATCH.EQ.0) THEN
      WRITE (*,*) 'THE FOLLOWING FILE ALREADY EXISTS:'
      WRITE (*,'(1X,A)') FNAME
      WRITE (*,*) 'ENTER THE NAME OF A NEW FILE (<CR>=QUIT):'
      READ (*,'(A)') FNAME
      IF (FNAME.EQ.' ') STOP
      GO TO 30
      ELSE
      IF (IO.EQ.1) WRITE (IOUT,*) 'FILE ALREADY EXISTS:'
      IF (IO.EQ.1) WRITE (IOUT,'(A)') FNAME
		call stopfile  ! emrl
      STOP
      END IF
      END IF
      FMT='FORMATTED'
      ACS='SEQUENTIAL'
      IF (IBIN.EQ.1) THEN
         FMT=FORM
         ACS=ACCESS
      END IF
      OPEN (IU,FILE=FNAME,STATUS='NEW',FORM=FMT,ACCESS=ACS,IOSTAT=IERR)
      IF (IERR.GT.0) GO TO 40
      RETURN
      END IF
C
C  OPEN A SCRATCH FILE
C
      IF (NSTAT.EQ.3) THEN
      FMT='FORMATTED'
      ACS='SEQUENTIAL'
      IF (IBIN.EQ.1) THEN
         FMT=FORM
         ACS=ACCESS
      END IF
      OPEN (IU,STATUS='SCRATCH',FORM=FMT,ACCESS=ACS,IOSTAT=IERR)
C  FOR MICROSOFT FORTRAN USE:
C      OPEN (IU,FORM=FMT,IOSTAT=IERR)
      IF (IERR.GT.0) GO TO 40
      RETURN
      END IF
C
C  OPEN A FILE OF UNKNOWN STATUS
C
      IF (NSTAT.EQ.4) THEN
      INQUIRE (FILE=FNAME,EXIST=EX,UNFORMATTED=FM)
      IF (EX) THEN
      FMT='FORMATTED'
      ACS='SEQUENTIAL'
      IF (IBIN.EQ.1.OR.FM.EQ.'YES') THEN
         FMT=FORM
         ACS=ACCESS
      END IF
      OPEN (IU,FILE=FNAME,STATUS='OLD',FORM=FMT,ACCESS=ACS,IOSTAT=IERR)
      IF (IERR.GT.0) GO TO 40
      ELSE
      FMT='FORMATTED'
      ACS='SEQUENTIAL'
      IF (IBIN.EQ.1) THEN
         FMT=FORM
         ACS=ACCESS
      END IF
      OPEN (IU,FILE=FNAME,STATUS='NEW',FORM=FMT,ACCESS=ACS,IOSTAT=IERR)
      IF (IERR.GT.0) GO TO 40
      END IF
      RETURN
      END IF
C
C  WRITE MESSAGE INDICATING PROBLEM OPENING FILE
C
40    IF (IBATCH.EQ.0) WRITE (*,5000) IU
      IF (IO.EQ.1 .AND. IOUT.GT.0) WRITE (IOUT,5000) IU
5000  FORMAT (' ERROR OPENING FILE TO UNIT ',I3)
	call stopfile  ! emrl
      STOP
      END
 
C***** SUBROUTINE *****
      SUBROUTINE CHOP(STRING,LNG)
      CHARACTER*(*) STRING
C
      LNGSTR= LEN(STRING)
      DO 1 N=LNGSTR,1, -1
      IF(STRING(N:N).NE.' ') THEN
      LNG=N
      GO TO 5
      END IF
1     CONTINUE
      LNG=0
5     RETURN
      END
 
C***** SUBROUTINE *****
      SUBROUTINE UPCASE(STRING)
      CHARACTER*(*) STRING
      CHARACTER*26 LETTER
      CHARACTER*1 CLC
C
      LETTER= 'ABCDEFGHIJKLMNOPQRSTUVWXYZ'
C
      IOFF= ICHAR('a') - ICHAR('A')
      LNG= LEN(STRING)
C
      DO 10 N=1,LNG
      DO 1 L=1,26
      ICLC= ICHAR(LETTER(L:L)) + IOFF
      CLC= CHAR(ICLC)
      IF(STRING(N:N).EQ.CLC) THEN
      STRING(N:N)=LETTER(L:L)
      GO TO 10
      END IF
1     CONTINUE
10    CONTINUE
C
      RETURN
      END
 
C***** SUBROUTINE *****
      SUBROUTINE UREADI(IARRAY,NDIM,IU,I7)
      DIMENSION IARRAY(NDIM)
      CHARACTER*133 LINE
C
      N1=1
5     READ(IU,'(A)') LINE
      ICOL=1
      DO 10 N=N1,NDIM
      CALL URWORD(LINE,ICOL,IWSTRT,IWLAST,2,IARRAY(N),RDUMMY,I7,IU)
      IF(LINE(IWSTRT:IWLAST).EQ.' ') THEN
      N1=N
      GO TO 5
      END IF
10    CONTINUE
C
      RETURN
      END
 
C***** SUBROUTINE *****
      SUBROUTINE ZERO(X,NJ,NI,NK)
      DIMENSION X(NJ,NI,NK)
      DO 10 K=1,NK
      DO 10 I=1,NI
      DO 10 J=1,NJ
10    X(J,I,K)=0.0
      RETURN
      END
 
C***** SUBROUTINE *****
      SUBROUTINE BALNCE (QX,QY,QZ,QSS,QSTO,IBOUND,HEAD,NCOL,NROW,NLAY,
     1NCP1,NRP1,NLP1,KOUNT,I7,ISILNT,ETOL,KPER,KSTP,ISS,HDRY)
C
      DIMENSION QX(NCP1,NROW,NLAY),QY(NCOL,NRP1,NLAY),
     1QZ(NCOL,NROW,NLP1),QSS(NCOL,NROW,NLAY),HEAD(NCOL,NROW,NLAY),
     2IBOUND(NCOL,NROW,NLAY),QSTO(NCOL,NROW,NLAY)
      DIMENSION IERRS(500),JERRS(500),KERRS(500),ERRS(500),ERRSA(500)
C
      NERRMX=500
      KOUNT=0
      EMAX=0.0
      DO 10  K=1,NLAY
      DO 10  I=1,NROW
      DO 10  J=1,NCOL
      IF(IBOUND(J,I,K).LE.0) GO TO 10
      IF(HEAD(J,I,K).EQ.HDRY) GO TO 10
      QX1=QX(J,I,K)
      QX2=QX(J+1,I,K)
      QY1=QY(J,I+1,K)
      QY2=QY(J,I,K)
      QZ1=QZ(J,I,K+1)
      QZ2=QZ(J,I,K)
C
      QI=0.0
      QO=0.0
      IF(QX1.GT.0.0) THEN
      QI=QI+QX1
      ELSE
      QO=QO-QX1
      END IF
      IF(QX2.LT.0.0) THEN
      QI=QI-QX2
      ELSE
      QO=QO+QX2
      END IF
      IF(QY1.GT.0.0) THEN
      QI=QI+QY1
      ELSE
      QO=QO-QY1
      END IF
      IF(QY2.LT.0.0) THEN
      QI=QI-QY2
      ELSE
      QO=QO+QY2
      END IF
      IF(QZ1.GT.0.0) THEN
      QI=QI+QZ1
      ELSE
      QO=QO-QZ1
      END IF
      IF(QZ2.LT.0.0)  THEN
      QI=QI-QZ2
      ELSE
      QO=QO+QZ2
      END IF
      IF(QSS(J,I,K).GT.0.0) QI=QI + QSS(J,I,K)
      IF(QSS(J,I,K).LT.0.0) QO=QO - QSS(J,I,K)
      IF(QSTO(J,I,K).GT.0.0) QI=QI + QSTO(J,I,K)
      IF(QSTO(J,I,K).LT.0.0) QO=QO - QSTO(J,I,K)
      QIMO= QI-QO
      QAVE=(QI+QO)/2.0
      IF(QAVE.GT.1.0E-20) THEN
      ERROR= 100.*(QIMO/QAVE)
      ELSE
      GO TO 10
      END IF
      ABERR= ABS(ERROR)
      ABEMAX= ABS(EMAX)
      IF(ABERR.GT.ABEMAX) THEN
        EMAX=ERROR
        IEMAX=I
        JEMAX=J
        KEMAX=K
      END IF
      IF(ABERR.GT.ETOL) THEN
        KOUNT=KOUNT+1
        IF(KOUNT.LE.NERRMX) THEN
          IERRS(KOUNT)=I
          JERRS(KOUNT)=J
          KERRS(KOUNT)=K
          ERRS(KOUNT)=ERROR
          ERRSA(KOUNT)=QIMO
        END IF
      END IF
10    CONTINUE
      IF(ISS.EQ.0) THEN
      IF (ISILNT.EQ.0) WRITE (*,5020) KPER,KSTP,KOUNT,ETOL,EMAX,IEMAX,
     1                                JEMAX,KEMAX
      WRITE(I7,*) ' '
      WRITE (I7,5020) KPER,KSTP,KOUNT,ETOL,EMAX,IEMAX,JEMAX,KEMAX
      ELSE IF(ISS.EQ.1) THEN
      IF (ISILNT.EQ.0) WRITE (*,5021) KOUNT,ETOL,EMAX,IEMAX,JEMAX,KEMAX
      WRITE(I7,*) ' '
      WRITE (I7,5021) KOUNT,ETOL,EMAX,IEMAX,JEMAX,KEMAX
      END IF
 
C  WRITE CELL-BY-CELL ERRORS THAT EXCEEDED TOLERANCE
      IF(KOUNT.GT.NERRMX) THEN
        NN=NERRMX
      ELSE
        NN=KOUNT
      END IF
      IF(NN.GT.0) THEN
        WRITE(I7,*)
     1  ' CELL-BY-CELL ERROR SUMMARY:  ROW  COLUMN  LAYER  ERROR(%)  ERR
     2OR(ABSOLUTE)'
        IF(KOUNT.GT.NERRMX)
     1    WRITE(I7,*) ' (ONLY FIRST 500 ERRORS ARE LISTED)'
        DO 15 N=1,NN
        WRITE(I7,'(3I5,F9.4,E15.5)') IERRS(N),JERRS(N),KERRS(N),ERRS(N),
     1                         ERRSA(N)
15      CONTINUE
      END IF
 
      RETURN
5020  FORMAT(1X,'PERIOD',I3,', STEP',I3,':',I6,' CELLS HAD ERRORS >',
     1F9.4,'%'/3X,'MAXIMUM ERROR = ',F9.4,'% IN ROW',I5,' COL',I5,
     2' LAYER',I5)
5021  FORMAT(1X,I6,' CELLS HAD ERRORS >',F9.4/
     1 3X,'MAXIMUM ERROR = ',F9.4,'% IN ROW',I5,' COL',I5,' LAYER',I5)
      END
 
C***** SUBROUTINE *****
      SUBROUTINE CELDAT (XMAX,YMAX,DELR,DELC,ZBOT,ZTOP,DELZCB,POR,
     1IBOUND,HEAD,QX,QY,QZ,LAYCON,NCOL,NROW,NLAY,NZDIM,IGRID,NLPOR,
     2NCP1,NRP1,NLP1,NCON,QSS,IBATCH,I7,QSTO,ISS,BUFF,NUMTS,
     3NPER,I10,IBSTRT)
C
      CHARACTER*78 MES
C
      DIMENSION IARRAY(10)
C
      DIMENSION XMAX(NCOL),YMAX(NROW),DELR(NCOL),DELC(NROW),ZBOT(NZDIM),
     1ZTOP(NZDIM),DELZCB(NLAY),POR(NCOL,NROW,NLPOR),
     2IBOUND(NCOL,NROW,NLAY),HEAD(NCOL,NROW,NLAY),QX(NCP1,NROW,NLAY),
     3QY(NCOL,NRP1,NLAY),QZ(NCOL,NROW,NLP1),NCON(NLAY),LAYCON(NLAY),
     4QSS(NCOL,NROW,NLAY),QSTO(NCOL,NROW,NLAY),NUMTS(NPER),
     5IBSTRT(NCOL,NROW,NLAY)
C
      NOLD=0
      MES= ' DO YOU WANT TO CHECK DATA FOR ANOTHER CELL ?'
      WRITE(I7,*) ' '
      WRITE(I7,*)
     1'***************************'
10    CONTINUE
      IF(ISS.EQ.0) THEN
      CALL MPRMPT
     1('ENTER: STRESS PERIOD & TIME STEP NUMBER')
      CALL GETI(2,IARRAY,0,0,0,'5.7.1')
      KPER=IARRAY(1)
      KSTP=IARRAY(2)
      CALL GETSTP(KPER,KSTP,NUMTS,NPER,I7,NSTEP)
      END IF
      CALL MPRMPT
     1('ENTER J I K COORDINATES OF CELL:  J   I   K')
      CALL GETI(3,IARRAY,0,0,0,'5.7.2')
      J=IARRAY(1)
      I=IARRAY(2)
      K=IARRAY(3)
      IF(ISS.EQ.0) THEN
      IF(NSTEP.NE.NOLD) THEN
      NOLD=NSTEP
      CALL READHQ(QX,QY,QZ,QSS,QSTO,HEAD,IBOUND,BUFF,NCOL,
     1 NROW,NLAY,NCP1,NRP1,NLP1,I10,NSTEP,DELT,TOTSIM,IERR,
     2 NUMTS,NPER,NSTEPF,NSTEPL,KPER,KSTP,1,I7)
      IF(IERR.NE.0) THEN
        IF(IBATCH.EQ.0) WRITE(*,5002) KPER,KSTP
        WRITE(I7,5002) KPER,KSTP
5002    FORMAT(1X,'PERIOD ',I4,' STEP ',I4,' NOT INCLUDED IN CBF.')
        GO TO 100
      END IF
      END IF
      END IF
      WRITE (I7,5000) I,J,K
      IF(ISS.EQ.0) WRITE (I7,5001) KPER,KSTP,DELT
      IF (IBATCH.EQ.0) WRITE (*,5000) I,J,K
      IF (IBATCH.EQ.0 .AND. ISS.EQ.0) WRITE (*,5001) KPER,KSTP,DELT
5000  FORMAT(' DATA FOR CELL: ROW',I4,', COLUMN',I4,', LAYER',I4)
5001  FORMAT(' STRESS PERIOD =',I4,', TIME STEP =',I4,
     1' TIME STEP LENGTH =',E13.5)
      IF(IBOUND(J,I,K).GT.-1000.AND.IBOUND(J,I,K).LT.1000) THEN
      IB=IBOUND(J,I,K)
      ELSE
      IB=IBOUND(J,I,K)/1000
      END IF
      IB= IB*IBSTRT(J,I,K)
      IF (IBATCH.EQ.0) WRITE (*,5010) IB,HEAD(J,I,K),POR(J,I,K)
      WRITE (I7,5010) IB,HEAD(J,I,K),POR(J,I,K)
5010  FORMAT(' IBOUND=',I7,' HEAD=',E15.5,' POROSITY=',F6.3)
      X1=0.0
      X2=XMAX(J)
      IF(J.GT.1) X1=XMAX(J-1)
      IF (IBATCH.EQ.0) WRITE (*,5020) X1,X2,DELR(J)
      WRITE (I7,5020) X1,X2,DELR(J)
5020  FORMAT(' X1=',E15.5,' X2=',E15.5,' (X2-X1)=',E15.5)
      Y1=0.0
      Y2=YMAX(I)
      IF(I.LT.NROW) Y1=YMAX(I+1)
      IF (IBATCH.EQ.0) WRITE (*,5030) Y1,Y2,DELC(I)
      WRITE (I7,5030) Y1,Y2,DELC(I)
5030  FORMAT(' Y1=',E15.5,' Y2=',E15.5,' (Y2-Y1)=',E15.5)
      IF(IGRID.EQ.1) THEN
      N=K
      ELSE
      N= (K-1)*NCOL*NROW + (I-1)*NCOL + J
      END IF
      ZMX=ZTOP(N)
      IF(LAYCON(K).EQ.1) ZMX=HEAD(J,I,K)
      IF(LAYCON(K).GT.1.AND.HEAD(J,I,K).LT.ZTOP(N)) ZMX=HEAD(J,I,K)
      DZ=ZMX-ZBOT(N)
      IF (IBATCH.EQ.0) WRITE (*,5040) ZBOT(N),ZMX,DZ
      WRITE (I7,5040) ZBOT(N),ZMX,DZ
5040  FORMAT(' Z1=',E15.5,' Z2=',E15.5,' (Z2-Z1)=',E15.5)
      QX1=QX(J,I,K)
      QX2=QX(J+1,I,K)
      IF (ABS(DZ).GT.1.0E-20) THEN
      VX1=QX1/(Y2-Y1)/(ZMX-ZBOT(N))/POR(J,I,K)
      VX2=QX2/(Y2-Y1)/(ZMX-ZBOT(N))/POR(J,I,K)
      ELSE
      VX1=0.0
      VX2=0.0
      END IF
      IF (IBATCH.EQ.0) THEN
      WRITE (*,*)
     1'VOLUMETRIC FLOW (Q1-Q6) AND VELOCITY (V1-V6) AT CELL FACES'
      WRITE (*,5050) QX1,QX2,VX1,VX2
5050  FORMAT(' X(FACES 1&2): Q1=',1PE13.5,' Q2=',E13.5,' V1=',E11.3,
     1 ' V2=',E11.3)
      END IF
      WRITE (I7,*)
     1'VOLUMETRIC FLOW (Q1-Q6) AND VELOCITY (V1-V6) AT CELL FACES'
      WRITE (I7,5050) QX1,QX2,VX1,VX2
      QY1=QY(J,I+1,K)
      QY2=QY(J,I,K)
      IF (ABS(DZ).GT.1.0E-20) THEN
      VY1=QY1/(X2-X1)/(ZMX-ZBOT(N))/POR(J,I,K)
      VY2=QY2/(X2-X1)/(ZMX-ZBOT(N))/POR(J,I,K)
      ELSE
      VY1=0.0
      VY2=0.0
      END IF
      IF (IBATCH.EQ.0) WRITE (*,5060) QY1,QY2,VY1,VY2
      WRITE (I7,5060) QY1,QY2,VY1,VY2
5060  FORMAT(' Y(FACES 3&4): Q3=',1PE13.5,' Q4=',E13.5,' V3=',E11.3,
     1 ' V4=',E11.3)
      IF(IB.NE.0) THEN
        QZ1=QZ(J,I,K+1)
        QZ2=QZ(J,I,K)
      ELSE
        QZ1=0.0
        QZ2=0.0
      END IF
      VZ1=QZ1/(X2-X1)/(Y2-Y1)/POR(J,I,K)
      VZ2=QZ2/(X2-X1)/(Y2-Y1)/POR(J,I,K)
      IF (IBATCH.EQ.0) WRITE (*,5070) QZ1,QZ2,VZ1,VZ2
      WRITE (I7,5070) QZ1,QZ2,VZ1,VZ2
5070  FORMAT(' Z(FACES 5&6): Q5=',1PE13.5,' Q6=',E13.5,' V5=',E11.3,
     1 ' V6=',E11.3)
      QI=0.0
      QO=0.0
      IF(QX1.GT.0.0) THEN
      QI=QI+QX1
      ELSE
      QO=QO-QX1
      END IF
      IF(QX2.LT.0.0) THEN
      QI=QI-QX2
      ELSE
      QO=QO+QX2
      END IF
      IF(QY1.GT.0.0) THEN
      QI=QI+QY1
      ELSE
      QO=QO-QY1
      END IF
      IF(QY2.LT.0.0) THEN
      QI=QI-QY2
      ELSE
      QO=QO+QY2
      END IF
      IF(QZ1.GT.0.0) THEN
      QI=QI+QZ1
      ELSE
      QO=QO-QZ1
      END IF
      IF(QZ2.LT.0.0)  THEN
      QI=QI-QZ2
      ELSE
      QO=QO+QZ2
      END IF
C
      QIQSS=0.0
      QOQSS=0.0
      QIQSTO=0.0
      QOQSTO=0.0
      IF(QSS(J,I,K).GT.0.0)  QIQSS=   QSS(J,I,K)
      IF(QSS(J,I,K).LT.0.0)  QOQSS=  -QSS(J,I,K)
      IF(QSTO(J,I,K).GT.0.0) QIQSTO=  QSTO(J,I,K)
      IF(QSTO(J,I,K).LT.0.0) QOQSTO= -QSTO(J,I,K)
      QITOT=QI + QIQSS + QIQSTO
      QOTOT=QO + QOQSS + QOQSTO
      DELTAQ=QITOT - QOTOT
 
      IF (IBATCH.EQ.0) THEN
      WRITE(*,5084)
      WRITE(*,5081) QI,QO
      WRITE(*,5082) QIQSS,QOQSS
      WRITE(*,5083) QIQSTO,QOQSTO
      WRITE(*,5086)
      WRITE(*,5087) QITOT,QOTOT
      WRITE(*,5088) DELTAQ
5081  FORMAT(1X,' FLOW ACROSS FACES = ',1PE13.5,' FLOW ACROSS FACES = ',
     1E13.5)
5082  FORMAT(1X,'       NET SOURCES = ',1PE13.5,'         NET SINKS = ',
     1E13.5)
5083  FORMAT(1X,'           STORAGE = ',1PE13.5,'           STORAGE = ',
     1E13.5)
5084  FORMAT(1X,'           INFLOW    ',13X,  '             OUTFLOW')
5085  FORMAT(1X,'           ------    ',13X,  '             -------')
5086  FORMAT(1X,' ------------------- ',13X,  ' --------------------- ')
5087  FORMAT(1X,'      TOTAL INFLOW = ',1PE13.5,'     TOTAL OUTFLOW = ',
     1E13.5)
5088  FORMAT(1X,'IN - OUT = ',1PE13.2)
      END IF
      WRITE(I7,5084)
      WRITE(I7,5085)
      WRITE(I7,5081) QI,QO
      WRITE(I7,5082) QIQSS,QOQSS
      WRITE(I7,5083) QIQSTO,QOQSTO
      WRITE(I7,5086)
      WRITE(I7,5087) QITOT,QOTOT
      WRITE(I7,5088) DELTAQ
 
      QAVE=(QITOT+QOTOT)/2.0
      IF(QAVE.GT.1.0E-20) THEN
      ERROR= 100.0*((QITOT-QOTOT)/QAVE)
      IF (IBATCH.EQ.0) then
      WRITE (*,5090) ERROR
      END IF
      WRITE (I7,5090) ERROR
      ELSE
      IF (IBATCH.EQ.0) THEN
      WRITE (*,*) 'NO FLOW INTO OR OUT OF THIS CELL'
      WRITE (*,*) '  NO MASS BALANCE COMPUTED'
      END IF
      WRITE (I7,*) 'NO FLOW INTO OR OUT OF THIS CELL'
      WRITE (I7,*) '  NO MASS BALANCE COMPUTED'
      END IF
5090  FORMAT(' PERCENT ERROR =',F10.5)
      IF(NCON(K).GT.0) THEN
      LCON=NCON(K)
      NLP=NLAY+LCON
      VZCB=QZ1/(X2-X1)/(Y2-Y1)/POR(J,I,NLP)
      IF (IGRID.EQ.1) THEN
      THICK=DELZCB(K)
      ELSE IF (IGRID.EQ.0.AND.K.LT.NLAY) THEN
      N= (K-1)*NCOL*NROW + (I-1)*NCOL + J
      NBELOW= N + (NCOL*NROW)
      THICK= ZBOT(N) - ZTOP(NBELOW)
      END IF
      IF(IBATCH.EQ.0) WRITE (*,5100) THICK,POR(J,I,NLP),VZCB
      WRITE (I7,5100) THICK,POR(J,I,NLP),VZCB
5100  FORMAT(' CONFINING BED:'/
     1'  THICKNESS =',1PE13.5,'  POROSITY =',E13.5,'  VZ =',E13.5)
      END IF
C
      WRITE(I7,*) ' '
      WRITE(I7,*)
     1'***************************'
100   CONTINUE
      CALL YESNO(MES,IANS,'5.7.3')
      IF(IANS.EQ.1) GO TO 10
      RETURN
      END
 
C***** SUBROUTINE *****
      SUBROUTINE FLIPQ(QX,QY,QZ,NCOL,NROW,NLAY,NCP1,NRP1,NLP1)
      DIMENSION QX(NCP1,NROW,NLAY),QY(NCOL,NRP1,NLAY),
     1 QZ(NCOL,NROW,NLP1)
C
C  SWITCH SIGNS OF FLOWS IF REVERSE TRACKING IS USED
C
      DO 10 K=1,NLAY
      DO 10 I=1,NROW
      DO 10 J=1,NCP1
10    QX(J,I,K)= -QX(J,I,K)
      DO 20 K=1,NLAY
      DO 20 I=1,NRP1
      DO 20 J=1,NCOL
20    QY(J,I,K)= - QY(J,I,K)
      DO 30 K=1,NLP1
      DO 30 I=1,NROW
      DO 30 J=1,NCOL
30    QZ(J,I,K)= -QZ(J,I,K)
      RETURN
      END
 
C***** SUBROUTINE *****
      SUBROUTINE GETPS(PERLEN,NUMTS,TIMX,NPER,TIMREL,TBEGIN,KKPER,
     1                 KKSTP,IERR)
      DIMENSION PERLEN(NPER),NUMTS(NPER),TIMX(NPER)
C
      IERR=0
      T=TBEGIN
      IF(TIMREL.LT.T) THEN
      IERR=1
      RETURN
      END IF
      DO 10 KP=1,NPER
      NSTPS=NUMTS(KP)
      TMULT=TIMX(KP)
      PERL=PERLEN(KP)
      DO 10 KS=1,NSTPS
      CALL STPSIZ(PERL,KS,NSTPS,TMULT,DT)
      TOLD=T
      T=T+DT
      IF(TIMREL.LE.T) THEN
      KKPER=KP
      KKSTP=KS
      TIMREL= (TIMREL-TOLD)/DT
      IF(TIMREL.LT.0.0) TIMREL=0.0
      IF(TIMREL.GT.1.0) TIMREL=1.0
      RETURN
      END IF
10    CONTINUE
      IERR=2
      RETURN
      END
 
C***** SUBROUTINE *****
      SUBROUTINE GETREQ(ITMS,ITREAD,MAXSTP,KOUNT,DELT,TFAC,T,I8)
C
C  ITMS = 0 -- LOCATIONS ARE NOT COMPUTED AT SPECIFIC POINTS IN TIME
C  ITMS = 1 -- LOCATIONS ARE COMPUTED AT SPECIFIC POINTS IN TIME
C
      IF (ITMS.EQ.1) THEN
C
C  ITREAD = 0 -- CONSTANT TIME STEP LENGTH
C  ITREAD = 1 -- TIME STEP DATA READ FROM A FILE
C
      IF (ITREAD.EQ.0) THEN
        IF(KOUNT.LE.MAXSTP) THEN
        T=T+DELT
        ELSE
        T=1.0E+29
        ITMS=0
        END IF
      ELSE IF (ITREAD.EQ.1) THEN
        IF(KOUNT.LE.MAXSTP) THEN
        READ(I8,*) T
        T= T*TFAC
        ELSE
        T=1.0E+29
        ITMS=0
        END IF
      END IF
 
      ELSE
      T= 1.0E+29
      END IF
      RETURN
      END
 
C***** SUBROUTINE *****
      SUBROUTINE GETSTP(KPER,KSTP,NUMTS,NPER,IU,NSTEP)
      DIMENSION NUMTS(NPER)
C
      NSTEP=0
      IF(KPER.GT.NPER) GO TO 20
      DO 10 K=1,KPER
      IF(K.EQ.KPER) THEN
      NSTP=KSTP
      ELSE
      NSTP=NUMTS(K)
      END IF
      IF(NSTP.GT.NUMTS(K)) GO TO 20
      DO 10 N=1,NSTP
      NSTEP=NSTEP+1
10    CONTINUE
      RETURN
20    WRITE(*,*)
     1'ERROR COMPUTING TIME STEP NUMBER IN ROUTINE <GETSTP>. STOP.'
      WRITE(IU,*)
     1'ERROR COMPUTING TIME STEP NUMBER IN ROUTINE <GETSTP>. STOP.'
		call stopfile  ! emrl
      STOP
      END
 
C***** SUBROUTINE *****
      SUBROUTINE STPSIZ(PERLEN,KSTP,NSTP,TMULT,STPL)
C
      IF(TMULT.NE.1.) THEN
      STPL=PERLEN*(1.-TMULT)/(1.-TMULT**NSTP)
      IF (KSTP.GT.1) THEN
      DO 10 N=2,KSTP
      STPL=TMULT*STPL
10    CONTINUE
      END IF
      ELSE
      RSTP=NSTP
      STPL=PERLEN/RSTP
      END IF
C
      RETURN
      END
 
C***** SUBROUTINE *****
      SUBROUTINE GENPAR(N,NRGEN,TBGABS,TBEGIN,TIMINI,TIMFIN,DTR,
     1  ICHEAD,ILOC,IBATCH,ILSTOR,NI,NJ,NK,IF1,IF2,IF3,IF4,IF5,IF6,
     2  NI1,NK1,NI2,NK2,NJ3,NK3,NJ4,NK4,NJ5,NI5,NJ6,NI6,NPART,NCOL,
     3  NROW,NLAY,ISS,KWT,IBOUND,HEAD,JLC,ILC,KLC,XLC,YLC,ZLLC,HNOFLO,
     4  HDRY,PERLEN,NUMTS,TIMX,TOT,INILOC,XMAX,YMAX,NPER,IUTEMP,
     5  IUSLOC,IUCBF,IUSUM,JMIN,JMAX,IMIN,IMAX,KMIN,KMAX,TSTOP)
      DIMENSION IBOUND(NCOL,NROW,NLAY),HEAD(NCOL,NROW,NLAY),
     1  JLC(NPART),ILC(NPART),KLC(NPART),XLC(NPART),YLC(NPART),
     2  ZLLC(NPART),TOT(NPART),PERLEN(NPER),NUMTS(NPER),TIMX(NPER),
     3  INILOC(NPART),XMAX(NCOL),YMAX(NROW)
 
      IPER=0
      ISTP=0
      IZRO=0
      DO 200 NR=1,NRGEN+1
C... CALCULATE VALUE OF THE REGENERATION TIME FOR THIS PASS THROUGH THE LOOP
      TRLEAS=TIMINI + (NR-1)*DTR
      IF(TRLEAS.GT.TIMFIN .AND. TIMFIN.GE.TIMINI) TRLEAS=TIMFIN
 
C... IF RELEASE TIME IS >= THE USER-SPECIFIED STOP TIME FOR THE TRACKING
C    SIMULATION, DO NOT GENERATE THE PARTICLES
      IF(TRLEAS.GE.TSTOP) GO TO 200
 
C... SET IDCODE TO 0 IF PARTICLE IS RELEASED AND ACTIVE AT TIME=0. IF IT
C    IS RELEASED LATER, SET IDCODE= -2
      IF(TRLEAS.EQ.0.0) THEN
        IDCODE=0
      ELSE
        IDCODE= -2
      END IF
 
C... IF RUN IS TRANSIENT, FIND OUT IF WE NEED TO READ HEADS SO THAT THEY ARE
C    CURRENT WITH THE REGENERATION TIME FOR THIS PASS THROUGH THE LOOP.
      IF(ISS.EQ.0) THEN
        T= TBGABS + TRLEAS
        CALL GETPS(PERLEN,NUMTS,TIMX,NPER,T,TBEGIN,KPER,KSTP,IERR)
        IF(IERR.GT.0) GO TO 200
        IF(KPER.NE.IPER .OR. KSTP.NE.ISTP) THEN
          IPER=KPER
          ISTP=KSTP
          CALL CBFHED(HEAD,NCOL,NROW,NLAY,IUCBF,NSTEP,IERR,NUMTS,NPER,
     1                NSTEPF,NSTEPL,KPER,KSTP,2,IUSUM)
        END IF
      END IF
 
      NOLD=N
      IF(NR.EQ.1) N1=NOLD
      DO 100 K=KMIN,KMAX
      DO 100 I=IMIN,IMAX
      DO 100 J=JMIN,JMAX
      KKK=K
 
C... IF KWT=1, LOOP THROUGH LAYERS TO FIND THE FIRST ACTIVE CELL THAT IS
C    ACTIVE AND WET.
      IF(KWT.EQ.1) THEN
      DO 110 KK=1,NLAY
      KKK=KK
      IF(HEAD(J,I,KKK).NE.HDRY .AND. HEAD(J,I,KKK).NE.HNOFLO) GO TO 120
110   CONTINUE
120   CONTINUE
      END IF
 
      ISKIP=0
C... SKIP OVER CELL IF IT IS DRY
      IF(HEAD(J,I,KKK).EQ.HNOFLO .OR. HEAD(J,I,KKK).EQ.HDRY) ISKIP=1
C... SKIP OVER CELL IF IT IS CONSTANT HEAD AND OPTION IS SET TO SKIP
C    CONSTANT HEAD CELLS
      IF (ICHEAD.EQ.0.AND.IBOUND(J,I,KKK).LT.0) ISKIP=1
C... CHECK ISKIP FLAG AND SKIP CELL IF NECESSARY
      IF(ISKIP.EQ.1) GO TO 100
 
C... IF THE RUN IS TRANSIENT (ISS=0), GENERATE COMPUTE NEW ARRAYS OF
C       PARTICLES EACH PASS THROUGH THE REGENERATION LOOP.
C    IF RUN IS STEADY STATE (ISS=1), ONLY GENERATE NEW PARTICLES THE
C       FIRST PASS THROUGH THE REGENERATION LOOP. ON THE NEXT PASSES,
C       JUST REUSE THE ARRAYS GENERATED ON THE FIRST PASS.
 
C... IF RUN IS TRANSIENT OR THIS IS THE FIRST PASS ...
      IF(ISS.EQ.0 .OR. NR.EQ.1) THEN
C
C... IF ILOC=0, GENERATE 3-D ARRAY BY VOLUME
      IF (ILOC.EQ.0) THEN
      CALL VOLDIV(J,I,KKK,JLC,ILC,KLC,XLC,YLC,ZLLC,NI,NJ,NK,NPART,N,
     1IBATCH)
      ELSE
 
C... IF ILOC NOT = 0, GENERATE 2-D ARRAYS BY FACE
      IF(IF1.EQ.1) CALL FACDIV (J,I,KKK,JLC,ILC,KLC,XLC,YLC,ZLLC,1,
     1NI1,NK1,NPART,N,IBATCH)
      IF(IF2.EQ.1) CALL FACDIV (J,I,KKK,JLC,ILC,KLC,XLC,YLC,ZLLC,2,
     1NI2,NK2,NPART,N,IBATCH)
      IF(IF3.EQ.1) CALL FACDIV (J,I,KKK,JLC,ILC,KLC,XLC,YLC,ZLLC,3,
     1NJ3,NK3,NPART,N,IBATCH)
      IF(IF4.EQ.1) CALL FACDIV (J,I,KKK,JLC,ILC,KLC,XLC,YLC,ZLLC,4,
     1NJ4,NK4,NPART,N,IBATCH)
      IF(IF5.EQ.1) CALL FACDIV (J,I,KKK,JLC,ILC,KLC,XLC,YLC,ZLLC,5,
     1NJ5,NI5,NPART,N,IBATCH)
      IF(IF6.EQ.1) CALL FACDIV (J,I,KKK,JLC,ILC,KLC,XLC,YLC,ZLLC,6,
     1NJ6,NI6,NPART,N,IBATCH)
      END IF
 
      END IF
 
100   CONTINUE
 
C... IF THIS IS A STEADY STATE RUN WITH MULTIPLE RELEASE, JUST
C    DUPLICATE THE PARTICLE COORDINATES FROM THE FIRST PASS
      IF(ISS.EQ.1 .AND. NR.GT.1) THEN
        IF(N2.GT.N1) THEN
          DO 130 NN=N1+1,N2
          N=N+1
          JLC(N)=JLC(NN)
          ILC(N)=ILC(NN)
          KLC(N)=KLC(NN)
          XLC(N)=XLC(NN)
          YLC(N)=YLC(NN)
          ZLLC(N)=ZLLC(NN)
130       CONTINUE
        END IF
      END IF
 
      IF(NR.EQ.1) N2=N
      IF(N.GT.NOLD) THEN
      DO 150 NN=NOLD+1,N
      INILOC(NN)= (KLC(NN)-1)*NCOL*NROW + (ILC(NN)-1)*NCOL + JLC(NN)
      IF(IDCODE.EQ.-2) THEN
        TOT(NN)= -1.0E+30
      ELSE
        TOT(NN)= -1.0
      END IF
 
C... SAVE STARTING LOCATIONS TO A FILE IF THAT OPTION WAS SELECTED
      IF(ILSTOR.EQ.1) WRITE(IUSLOC,5210) JLC(NN),ILC(NN),KLC(NN),
     1                  XLC(NN),YLC(NN),ZLLC(NN),0,0,0,TRLEAS
 
C... CONVERT TO GLOBAL X & Y COORDINATES IF NEEDED
      IF(ISS.EQ.0 .OR. NR.EQ.1) THEN
        JP=JLC(NN)
        IP=ILC(NN)
        XMN=0.0
        XMX=XMAX(JP)
        IF(JP.GT.1) XMN= XMAX(JP-1)
        XLC(NN)= (1.0-XLC(NN))*XMN + XLC(NN)*XMX
        YMN=0.0
        YMX=YMAX(IP)
        IF(IP.LT.NROW) YMN=YMAX(IP+1)
        YLC(NN)= (1.0-YLC(NN))*YMN + YLC(NN)*YMX
      END IF
 
C... WRITE SOME PARTICLE LOCATION DATA TO A SCRATCH DIRECT ACCESS FILE
      WRITE(IUTEMP,REC=NN) XLC(NN),YLC(NN),ZLLC(NN),IDCODE,TRLEAS,IZRO
150   CONTINUE
      END IF
 
200   CONTINUE
      RETURN
 
5210  FORMAT(3I5,2X,3F12.7,2X,3(I1,1X),1X,E13.5)
      END
 
C***** SUBROUTINE *****
      SUBROUTINE GETABT(PERLEN,NUMTS,TIMX,NPER,TIMREL,TBEGIN,TIMAB,KPER,
     1                 KSTP,IERR)
      DIMENSION PERLEN(NPER),NUMTS(NPER),TIMX(NPER)
C
      IERR=0
      IF(TIMREL.LT.0.0 .OR. TIMREL.GT.1.0) IERR=1
      IF(KPER.GT.NPER) IERR=1
      IF(KSTP.GT.NUMTS(KPER)) IERR=1
      IF(IERR.NE.0) RETURN
 
      T=TBEGIN
      DO 10 KP=1,KPER
      KSPER=NUMTS(KP)
      NSTPS=KSPER
      IF(KP.EQ.KPER) NSTPS=KSTP
      TMULT=TIMX(KP)
      PERL=PERLEN(KP)
      DO 10 KS=1,NSTPS
      CALL STPSIZ(PERL,KS,KSPER,TMULT,DT)
      TOLD=T
      T=T+DT
10    CONTINUE
C
      TIMAB= TIMREL*T + (1.0-TIMREL)*TOLD
C
      RETURN
      END
 
C***** SUBROUTINE *****
      SUBROUTINE ND2IJK(ND,I,J,K,NROW,NCOL)
      NRC= NROW*NCOL
      K= 1 + (ND-1)/NRC
      I= 1 + (ND - (K-1)*NRC - 1)/NCOL
      J= ND - (K-1)*NRC - (I-1)*NCOL
      RETURN
      END
 
C***** SUBROUTINE *****
      SUBROUTINE RELEAS(NP,TOT,ZTOP,ZBOT,HEAD,LAYCON,ZLC,ZLLC,NDFRST,
     1  NZDIM,NCOL,NROW,NLAY,TIME,TMAX,IUPART,IUSUM,IGRID,HDRY,HNOFLO,
     2  IPSTAT,NSTEP)
      DIMENSION ZTOP(NZDIM),ZBOT(NZDIM),HEAD(NCOL,NROW,NLAY),
     1  LAYCON(NLAY)
C
      IPSTAT=0
C
      READ(IUPART,REC=NP) XSTRT,YSTRT,ZLSTRT,IDCODE,TRLEAS,NS
 
C... IF THE RELEASE TIME IS LESS THAN THE MINIMUM STARTING TIME FOR
C    THIS PASS THROUGH THE TIME LOOP, RETURN.
      IF(TRLEAS.LT.TIME) RETURN
C... IF THE RELEASE TIME IS >= TMAX, SET IPSTAT=1 AND RETURN.
      IF(TRLEAS.GE.TMAX) THEN
        IPSTAT=1
        RETURN
      END IF
 
C... IF PARTICLE IS ALREADY ACTIVE, SET TOT= -1.0 & RETURN
      IF(IDCODE.EQ.0) THEN
        TOT= -1.0
        RETURN
      END IF
 
      CALL ND2IJK(NDFRST,I,J,K,NROW,NCOL)
 
C... IF CELL IS DRY OR INACTIVE, JUST RETURN
      IF(HEAD(J,I,K).EQ.HDRY .OR. HEAD(J,I,K).EQ.HNOFLO) RETURN
 
      IF(IGRID.EQ.1) THEN
        TOP=ZTOP(K)
        BOT=ZBOT(K)
      ELSE
        TOP=ZTOP(NDFRST)
        BOT=ZBOT(NDFRST)
      END IF
 
      IF(LAYCON(K).EQ.1) TOP=HEAD(J,I,K)
      IF(LAYCON(K).GT.1 .AND. HEAD(J,I,K).LT.TOP) TOP=HEAD(J,I,K)
 
 
C... IF THE LOCAL Z COORDINATE IS VALID, THEN COMPUTE GLOBAL Z COORDINATE
C    IF NECESSARY. RELEASE THE PARTICLE IF THE GLOBAL COORDINATED IS
C    WITHIN THE SATURATED THICKNESS OF THE CELL. OTHERWISE, JUST RETURN.
C
      IF(ZLLC.GE.-1.0 .AND. ZLLC.LE.1.0) THEN
          CALL ZCALC(IGRID,NCOL,NROW,NLAY,Z,ZLLC,ZTOP,ZBOT,HEAD(J,I,K),
     1               LAYCON(K),I,J,K,NZDIM,IUSUM)
        IF(Z.GT.TOP .OR. Z.LT.BOT) RETURN
        ZLC=Z
        GO TO 100
C
C... ELSE IF THE LOCAL COORDINATE IS SET TO 1.0E+30, CALCULATE THE LOCAL
C    COORDINATE FROM THE GLOBAL COORDINATE. IF THE LOCAL COORDINATE IS
C    WITHIN THE SATURATED THICKNESS OF THE CELL (0 TO 1), RELEASE THE
C    PARTICLE. OTHERWISE, JUST RETURN.
C
      ELSE IF(ZLLC.EQ.1.0E+30) THEN
        IF(ZLC.EQ.1.0E+30) RETURN
        Z= (ZLC-BOT)/(TOP-BOT)
        IF(Z.GT.1.0 .OR. Z.LT.0.0) RETURN
        IF(Z.EQ.1.0) THEN
          Z=0.999
          ZLC= Z*TOP + (1.0-Z)*BOT
        ELSE IF(Z.EQ.0.0) THEN
          Z=0.001
          ZLC= Z*TOP + (1.0-Z)*BOT
        END IF
        ZLLC=Z
        GO TO 100
C
C... IF IT GETS THIS FAR, THERE IS SOMETHING WRONG WITH THE PARTICLE'S
C    COORDINATES, SO JUST RETURN.
C
      ELSE
        RETURN
      END IF
 
C... IF THE PARTICLE WAS ACTIVATED, UPDATE PARAMETERS AND REWRITE THEM
C    TO THE PARTICLE SCRATCH FILE
100   CONTINUE
      IDCODE=0
      TOT= -1.0
      ZLSTRT=ZLLC
      TIME=TRLEAS
      WRITE(IUPART,REC=NP) XSTRT,YSTRT,ZLSTRT,IDCODE,TRLEAS,NSTEP
      RETURN
      END
 
C***** SUBROUTINE *****
      SUBROUTINE ZCALC (IGRID,NCOL,NROW,NLAY,ZLC,ZLOC,ZTOP,ZBOT,
     1 HEAD,LAYCON,IP,JP,KP,NZDIM,IUSUM)
C
      DIMENSION ZTOP(NZDIM),ZBOT(NZDIM)
C
      IF(IGRID.EQ.1) THEN
      NK=KP
      NKP=NK+1
      ELSE
      NK = (KP-1)*NCOL*NROW + (IP-1)*NCOL + JP
      NKP= NK + (NCOL*NROW)
      END IF
 
      IF(KP.EQ.NLAY.AND.ZLOC.LT.0.0) THEN
      WRITE(IUSUM,*)' RUN STOPPED. INCONSISTENT CONFINING BED LOCATION'
		call stopfile  ! emrl
      STOP
      END IF
 
      ZMX=ZTOP(NK)
      IF(LAYCON.EQ.1) ZMX=HEAD
      IF(LAYCON.GT.1.AND.HEAD.LT.ZMX) ZMX=HEAD
      IF(ZLOC.GE.0.0) THEN
      ZLC=ZLOC*ZMX + (1.0-ZLOC)*ZBOT(NK)
      ELSE
      ZL= 1.0+ZLOC
      ZLC=ZL*ZBOT(NK) + (1.0-ZL)*ZTOP(NKP)
      END IF
      RETURN
      END
